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Complex networks exhibit a wide range of collective dynamic phenomena, including synchroniza¬ 
tion, diffusion, relaxation, and coordination processes. Their asymptotic dynamics is generically 
characterized by the local Jacobian, graph Laplacian or a similar linear operator. The structure of 
networks with regular, small-world and random connectivities are reasonably well understood, but 
their collective dynamical properties remain largely unknown. Here we present a two-stage mean- 
field theory to derive analytic expressions for network spectra. A single formula covers the spectrum 
from regular via small-world to strongly randomized topologies in Watts-Strogatz networks, explain¬ 
ing the simultaneous dependencies on network size N, average degree k and topological randomness 
q. We present simplified analytic predictions for the second largest and smallest eigenvalue, and 
numerical checks confirm our theoretical predictions for zero, small and moderate topological ran¬ 
domness q, including the entire small-world regime. For large q of the order of one, we apply 
standard random matrix theory thereby overarching the full range from regular to randomized net¬ 
work topologies. These results may contribute to our analytic and mechanistic understanding of 
collective relaxation phenomena of network dynamical systems. 

PACS numbers: 89.75.-k, 05.45.Xt, 87.19.1m 


I. INTRODUCTION 

The structural features of complex networks underlie 
their collective dynamics such as synchronization, diffu¬ 
sion, relaxation, and coordination processes Hi. Such 
processes occur in various fields, ranging from opinion 
formation in social networks i and consensus dynamics 
of agents [|j to synchronization in biological circuits 
and oscillations in gene regulatory networks and neural 
circuits The asymptotic collective dynamics of all 

such processes is characterized by the local Jacobian, the 
graph Laplacian or similar linear operators. In general, 
such linear operators directly connect the structure of an 
underlyin g ne twork to its dynamics and thus its function 
(see e.g. pjj, HU)- Therefore, a broad area of research 
is related to the study of properties of such operators, in 
particular to the study of their spectra [12Ml9l | . 

Small-world models based on rewiring have received 
widespread attention both theoretically and in applica¬ 
tions, as demonstrated, for instance, by the huge number 
of references pointing to the original theoretical work [20[ . 
But for most of their features analytical predictions are 
not known to date (fUJ; a mean held solution of its aver¬ 
age path length constitutes a notable exception [22j]). In 
particular, the spectrum of small-world Laplacians has 
been studied for several specific cases and numerically 
f23U27| yet a general derivation of reliable analytic pre¬ 


dictions was missing. 

Here we present a mean held theory [28| towards clos¬ 
ing this gap. The article is organized as follows. In Sec¬ 
tion nu we hrst review the relations between relaxation 
dynamics and the spectrum of the graph Laplacian. In 
Section fllll we present rewiring ‘on average’, a new mean 
held rewiring recently proposed in a brief report (28j . 
Based on this rewiring, we derive a single formula that 
approximates well the entire spectrum from regular to 
strongly randomized topologies. We then investigate the 
ordering of the mean held eigenspectrum in Section IIVI 
In Section [V] we quantify the accuracy of our predictions 
via systematic numerical checks for the extreme eigen¬ 
values. For the topological randomness q of the order of 
unity, standard random matrix theory is applied in Sec¬ 
tion IVII We close in Section IVIII with a summary and a 
discussion of further work. 


II. NETWORK RELAXATION DYNAMICS 

The relaxation dynamics towards equilibrium and re¬ 
lated collective phenomena close to invariant sets or sta¬ 
tionary distributions emerge across a wide variety of sys¬ 
tems ubiquitous in natural and artihcial systems [lj, |29j . 
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A. Generic linear relaxation 


B. Different example systems 


Mathematically, the relaxation of network dynamics 
to a fixed point, a periodic orbit or a similar stationary 
regime is generically characterized by equations of the 
form 


dxi 

dt 


N 

^ ", Jij ( x j ~ x i) 

j=l 


fori £ {1, ...,1V} , 


(HI) 


where Xi(t) quantifies the deviation at time t from an 
invariant state, IV £ N is the size of the network and 
Jij £ R quantifies the influence of unit j onto unit i. In 
general Xi(t) £ R d , we here take d = 1 for simplicity of 
presentation. 

The equivalent mathematical description 


dxi 

dt 


N 

Aij x j 

3 =1 


(n.2) 


for i £ {1, ...,1V} with the Laplacian 


Aij — 


Jij 


for i ^ j 


- Ef=i Jij for * = 3 


(II.3) 


follows directly from the original dynamics (III. II) . 

The eigenvalues A ra £ C and corresponding eigenvec¬ 
tors v n of such a Laplacian, satisfying 


Av n — A n v n (II.4) 

for n £ {1,..., IV}, fully characterize the asymptotic 
(linear or linearized) dynamics. For instance, for sta¬ 
ble dynamics, where all Xi(t) -A 0 for t -A oo, the 
largest nonzero (principal) eigenvalue A* dominates the 
long term dynamics: if we have 2 ^( 0 ) = y} n _ 1 a n v n , the 
vector x(t) = (xi(t),... ,Xat(1)) t evolves as 


x(t) = exp(Ai)x(0) (II.5) 

= exp (At) ^2 a nVn 

n 

N 

= ^ anexp(A n l)u n . 

n =1 

Due to stability we have a n = 0 whenever A„ = 0, and 
for long times this is dominated by 

x(t) ~ a* exp(A»t)n*. (II.6) 

where A* is the principal eigenvalue. 

Note that (III.5I) also reveals how exactly all eigenvalues 
contribute to relaxation (and how much relative to each 
other). Additional individual eigenvalues of interest are 
given by the one with smallest real part A_, because it 
bounds the real parts of the spectrum below and thereby 
determines the support of the spectrum, and also because 
it is involved in determining synchronizability conditions 
in coupled chaotic systems via the ratio A„/A_, see [30l| . 


We briefly consider two very different paradigmatic ex¬ 
ample classes of systems and comment on a few others 
whose relaxation properties are characterized by equa¬ 
tions of the same type as (III. 21) . 

Stochastic processes. Firstly, consider random walks 
on a graph, or equivalently, Markov chains defined by a 
weighted non-negative graph whose nodes represent the 
N states (3lj. For such processes, the dynamics of the 
probability Piit) to reside state i £ {1,..., N} at time t 
is given by 


dpi 

dt 


N 


yz i^pj - 7 /-A 

j= i 


fori £ {1, ...,1V} , (II.7) 


where defines the transition rate (probability per unit 
time) of the system switching state from j to i given 
it resides in i. We assume the process to be ergodic. 
Identifying p* to be the unique stationary distribution 
Ajj = Tij for i ^ j and An = — £W Tji an d setting 
Xi(t) = piit ) — p* exactly maps this processes onto the 
generic form (III.21) with Xi(t) —> 0 as t —> oo for all 
i £ {1,..., iV}. 

Coupled deterministic oscillators. Secondly, consider 
the relaxation dynamics of weakly coupled limit cycle os¬ 
cillators, generically modeled as phase-coupled oscillators 

+ for* £ {1,..., N} , (II.8) 
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where 0i(t) is the phase of unit i at time t, iOi is the 
local intrinsic frequency of oscillator i and hi A.) defines 
a smooth coupling function from unit j to i [y, 32, 33|. 
Phase-locking, where 6j(t) — 9i(t) = Aji is constant in 
time, constitutes a generic collective state of such sys¬ 
tems, see e.g. [34l - |36| . A paradigmatic model is given 
by networks of Kuramoto oscillators coupled by_simple 
sinusoidal functions, i.e. hij(9) = sin(0) |32l. l33i l37j. 

In the most general setting, a matrix J is defined by el¬ 
ements Jij = dhij (9)/d9\g—Aji that encode the network 
structure close to the phase-locked state. Under certain 
conditions on the w, and the hij, the system’s dynamics 
exhibits a short transient dominated by nonlinear effects 
and thereafter exponentially relaxes to the phase-locked 
state. Linearizing close to such a state yields phase per¬ 
turbations defined as 


Si(t) := 6i(t) - 9(t) (II.9) 

evolve according to 

f for i £ {1,..., IV} (11.10) 

j 

with the graph Laplacian given by (111.31) . 

In a simple setting, we have Jij = 1 for an existing 
edge and J^ = 0 for no edge such that the local linear 












operator in (III. 101) coincides with the graph Laplacian 
defined by its elements 
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A ij = J ij {l~8 ij )-k i 6 ij (11.11) 

for i,j £ {1, ...,1V}, where now Jij are the elements 
of the adjacency matrix, ki is the degree of node i (re¬ 
placed by the in-degree for directed networks) and <5y is 
the Kronecker-delta. The asymptotic relaxation dynam¬ 
ics on such networks is thus characterized by this graph 
Laplacian A. Similarly, any dynamics near genuine fixed 
points, for instance in gene regulatory networks {3§|, [39j 
is equally characterized by linearized dynamics stemming 
from local Jacobians. 

Power grids, social networks, neural circuits,... Sev¬ 
eral other systems exhibit qualitatively the same dynam¬ 
ical relaxation. In fact, power grids are often character¬ 
ized by second order oscillatory systems (40| that prin¬ 
cipally relax to periodic phase-locked solutions (station¬ 
ary operating states of the grid) very similarly to the 
phase oscillator systems discussed above [HI]. I n models 
of several social phenomena, e.g. of opinion formation of 
agents, the dynamics of consensus formation is equally 
akin to such locking dynamics, where the locked state 
would now be a homogeneous, fully synchronous one [42j. 
Last but not least, perturbations of the spatio-temporal 
collective dynamics of pulse-coupled systems such as neu¬ 
ral circuits [43, S3 also relax according to (III.2I) . 


/ w \ 

one random realization mean field 



FIG. 1. Rewiring — stochastic and mean field. (Cartoon 
for N = 10 and k = 4) Instead of taking out (step 1) and 
putting back (step 2) edges randomly (left column) the corre¬ 
sponding weight is subtracted (step 1) uniformly and added 
(step 2) in two fractions (right column). 


III. MEAN FIELD REWIRING AND 
SPECTRUM 

Diving into explaining the small-world model, we an¬ 
alyze and derive its approximate Laplacian based on a 
two stage mean field rewiring. We follow |[28j and where 
appropriate take parts of the description presented there. 

We consider an initial ring graph of N nodes. Each 
node receives k (being even) links from its k/2 nearest 
neighbors on both sides. Then we introduce randomness 
in the network topology by rewiring. 

To define Watts-Strogatz randomized networks, sin¬ 
gle instances of an ensemble of stochastically rewired 
networks are generated (Fig. [1] left panel). Following 
[20l ] for undirected networks, we first cut each edge with 
probability q. Afterwards the cut edges are rewired to 
nodes chosen uniformly at random from the whole net¬ 
work. Similarly, for directed [45| networks, we first cut 
all out-going edges with probability q and rewire their 
tips afterwards. In both cases we avoid double edges and 
self-loops. 

To analytically determine the Laplacian mean field 
spectrum in dependence of the network size N, the aver¬ 
age degree k and the topological randomness q, we intro¬ 
duce a two-stage mean field rewiring that effectively gen¬ 
erates, at given q the average network from the ensemble 
of all stochastically rewired networks. This is depicted in 
Fig- H (right panel) in comparison to both other rewiring 


procedures for undirected and directed networks. Firstly, 
we define a circulant mean field Laplacian 


( CO 

Cl 

C2 


Cat_i\ 

Cjv-l 

Co 

Cl 

C2 



Cat-1 

Co 

Cl 


^ Cl 




C 2 

Cl 

Cat-1 Co / 


Its matrix elements for the initial configuration (Fig. [T] 
q=0, top) are given by 

—k if i = 0 

0 if * e (| + 1,..., IV — | — 1} = 5 2 , 

(HI-2) 

where Si represents the set of edges being present in the 
ring and S 2 those absent ones outside that ring. 

Instead of rewiring single edges randomly (Fig. [T] 
left panel), we now distribute the corresponding weight 
of rewired edges uniformly among the whole network 
(Fig. |l] right panel). Thus, for a given rewiring probabil¬ 
ity q we generate a mean field version of the randomized 
network ensemble in the following two steps: 

Firstly, we subtract the average total weight qkN/2 of 
all edges to be rewired (Si), i.e. c, = 1 — q if i £ Si. 
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Secondly, the rewired weight is distributed uniformly 
among the total ’available’ weight in the whole network 
given by 


/ = «(«-D- (1-^JV . (m.3) 


With the weights 


h = 


qkN 


(in. 4) 


being available in Si and 

N(N - 1 )-kN 

h =- t, - 


(III- 5 ) 


in S' 2 , we assign the fraction f\/f to elements represent¬ 
ing edges in Si and / 2 // to those representing 5' 2 . There¬ 
fore, an individual edge in S\ gets the additional weight 


_ fi 2 _ g 2 ^ 

f hf N -l-(l-q)k' 


(HI 6) 



{ 

k/2 

}{N- 

-k 

- 1 }{ k/2 

} 

-k 

Wl • 


Wl W2 


W 2 Wl ■ ■ ■ ■ 

Wl 

Wl 

—k 

Wi - 

- ■ ■ Wl 

W 2 

♦ ♦ 

■ 

■ 

Wl 

—k 

Wl • ■ ■ 

Wl 

* ♦ 

• * ♦ 

Wl 

Wl 



* * * * 


*♦ * * 

W2 

W2 

♦, 



* ♦ , 

» * * # ♦ , 

. 


* # 

* * 



* ♦ * * 

W 2 

W2 

Wi 



* . Wl 


Wl —k uii 

Wl 

• 


* ♦ 

W2 

Wl 

■ ■ ■ ■ wi —k 

Wi 

Wl 

. . . . 

Wl 

W2 • ■ ■ 

1 W2 

Wl ■ ■ ■ ■ Wl 

-k 


FIG. 2. The banded structure of the mean field graph 
Laplacian A mt given in Eqs. (IIII.ll) and (IIII.8I) . It has 

the weights u)i = 1 — g + un for a\i £ Si and W 2 for a\i £ S 2 
(see Eq. (1III.2I) for the definition of Si). For q = 0 and hence 
wi = 1 and W 2 — 0 we recover the exact ring Laplacian. 


and an edge in S 2 gets the new weight 


VJ 2 = 


h 


qkN 

2 


qk 


f N(N-i)-kN JV — 1 — (1 — q)k 


(III.7) 


Thus, as depicted in Fig. [2l in our mean field theory the 
elements of the Laplacian d mf (IIII.ll) of a network on N 
nodes with degree k after rewiring with probability q are 
given by 


—k if i = 0 

Ci= { 1 — q + W\ if i £ Si 
W 2 if i £ S 2 ■ 


(III.8) 


The mean field Laplacian defined by (IIII.ll) and (IIII.8I) 
by construction is a circulant matrix with eigenvalues 


\ mf 
A l 


N—l 

C 1 exp 
1=0 


^ —27ri(Z - l)j ^ 


(III.9) 


Observing the structure in Fig. [5] we immediately ob¬ 
tain the trivial eigenvalue for 1 = 1: 


N -1 

A™ f = ^ Cj = —k + k( 1 — q + w 1 ) + (TV — k — 1 )w 2 = 0 , 
j=o 

(III.10) 

which is common to all networks (for all q , N, and any 
k < N — 1) and reflects the invariance of Laplacian dy¬ 
namics against uniform shifts, as seen from the associated 
eigenvector v\ = (1,..., 1) T . 

To obtain the remaining eigenvalues for l £ {2,..., IV} 
we first define 


Xi := exp 


1) ^ 


(III.ll) 


and 


c' := ^ + qc" (III.12) 


c" := q 


N — 1 — (1 — q)k 


(III.13) 


This leads to 


= - k + kc'^xj 
3 =1 

N~ 1 -| N—l 

+ kc" xj + kc' ^ xj 

1=4+1 3 = n -4 

k_ k_ 

= — k + kc' ^ x\ + kc' x f--* 
l=i 1=1 


(III. 14) 


+ fc" y . 


1=1 


Z x i 

1=1 


N_. N_ 

2 J +X, 2 


(III.15) 


where we have exploited the additional transposition 
symmetry A ml = which implies Cj = Cn-j ■ Ap¬ 

plying the Euler formula exp (ia) = cos ( a ) + i sin (a) , 
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FIG. 3. Interpolating the eigenvalues. Eqs. (IIII.18I) 

(blue) and (IIII.221) (red) both contain the eigenspectrum A™ f 
for l G {2,..., 2V} correctly (green circles). While A“ lf in 
Eq. (IIII. 1811 includes A) nf as lim/^i A) nf = A™ f = 0 as well, A™ f 
for l G {2,..., IV} in Eq. (IIII.2211 does not: To further simplify 
expressions, we have used identities (IIII. 191) and (IIII.201) only 
valid for integer l, but apparently not for l = 0. 


the complex summands cancel and we get 


k 



(III.16) 


Using the identity 


Y cos (ja) = 
j =o 


cos 


n + 1 


a sin 


(na\ 1 


V 2 7 sin(f) 


+ 1 


_ 1 (-. sin((n + \)a) 
2 l + sin(f) 


(III.17) 


we obtain 


sm 


A; = - k + kc' 


+ x, 2 kc"- 


( (fe+l)(Z —l)ir 

V w 


Sill (^l) 




in ( (AT-fc-iM-tU ) 

sin (Jtgz'j 


(III.18) 


Taking advantage of additional identities - only valid 
for l G Z (Fig. 0 


= (- 1 ) 1 - 1 , (ni.i9) 


(—1) ( 1 sin(a) = sin(a + (l — l)7r) (III.20) 

and the symmetry sin(— a) = — sin(a), the expression 
simplifies to 


\ mf 
A l 


— k + kc' 


— k + kc' 


— k + kc' 



- 1 I + (-ly-'fcc" 


sm 


(—(fc+l)+iV)(Z —1)TT 
N 


sm 


- 1 


kc 1. 


sin 

/ (— (fc+1) (Z — 1) +2JV(Z—1)) 7r 


N 


-1 - 


sin 

sm ( 


kc" ^ 


sm 




which finally leads to 


\mf 

A l 


—k — kc' + k (c ' 


c") 


sin ( (fc+iKi-iK ) 

(^ 


for / G {2,..., N}. 


(III.21) 


(III.22) 
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FIG. 4. A"' 1 always constitutes the second largest 

eigenvalue. Functions f(x) (IIV.3II . the oscillating function 
sin ((k + 1)*) and the envelope function 1/sin(x) are plotted 
vs. x = G (0,7r) for k = 10. Obviously, a larger k 

leads to more roots of f(x), but otherwise the functions show 
the same characteristics for all k < N — 1: f(x) has a local 
maximum at x = 0 and decreases strictly monotonically up 
to the following minimum. For larger x the envelope function 
guarantees that all values up to x = 7r/2 are smaller than 

f(xi =2 = #)■ 


IV. THE ORDERING OF THE MEAN FIELD 
SPECTRUM 


FIG. 5. Important points of the function f(x). The 
function f(x) is plotted for k = 100. The boundary points 
(green, eq. (IIV.9I) ') of the function f(x) and the envelope 
function 1/sin(a;), roots of f(x) (purple, eq. (IIV.8II 1 and the 
x-value Xi =2 (blue, eq. (1IV.7I) . N = 1000) corresponding to 
the eigenvalue A™* are highlighted. 


damped oscillation with period of 2ir/(k + 1) and with 
the amplitude decreasing as l/sin(x) (Fig. [4]). 

At x — 0 we apply the Theorem of l’Hospital to calcu¬ 
late the following limits. There is a removable singularity 

lim f( x ) = k + 1 (IV. 5) 

£—>■0 


The spectrum obeys the symmetry 

\f = \$_ l+2 , (iv. 1) 

but is unordered otherwise, i.e. the index l does neither 
denote eigenvalues with decreasing real part nor eigen¬ 
values with decreasing absolute value. 

As we argue below the expression A‘ 2 nf (that equals A™ f 
due to IlIV.ll) ) always constitutes the second largest (prin¬ 
cipal) eigenvalue A™ f . The only term depending on l in 
eq. (IIII.221) is the ratio 


with 

lim f'(x) = 0 and lim f"(x) = — \k{k 2 + 3k + 2) < 0 , 

(IV.6) 

i.e. a local maximum. 

In order to show that the index l = 2 is always asso¬ 
ciated with the second largest eigenvalue, we first deter¬ 
mine its x-value. It is given by 

*1=2 = £ ■ ( IVJ ) 


sin 
sin ( 


(IV.2) 


Since the roots of the function f[x) are located at 

7*7r 

•Eroot ,r ^ ^ (IV.8) 


We therefore study the function 


with 


fM _sm((k+» x ) _ 

(IV.3) 

sinx 

X - (* - 1 > 

x ~ N 

(IV.4) 


and x G (0,7t/2). Due to the symmetry (IIV.1I) the inter¬ 
val (0,7t/2) covers the entire spectrum (IIII.221) . 

The function f(x) on x G (0, 7t/2) is the product of the 
oscillating function sin((fc + l)x) and a strictly mono¬ 
tonically decreasing function 1/sin(x). Therefore, it is a 


for r G Z. Thus, Xi- 2 is always smaller than the first 
root Xroot.i (IIV.8I) of the function f(x) (Fig. [5]). 

The boundary points of function f(x) and the envelope 
function 1 / sin(x) are given by 


Anr + 7r 

Xb ’ r = 2(fcTT) ' 


(IV-9) 


for r G Z. 

The function f(x) is bounded from above by the enve¬ 
lope function 1 / sin(x) for all x > Xb,o with 


Xb, 0 


7r 

2 {k + 1) 


(IV.10) 
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being the first boundary point of function f(x ) and its 
envelope function (r = 0 in eq. (HV.9jl l (Fig. [5]). 

The first derivative of /( x) stays negative at least up 
to the first root (r = 1 in eq. (IIV. 81) 1 at 

7r 

•^root.l = ^ ^ ^ *^6,0 (IV.11) 

which is always larger than the first boundary point ccf, o 

divm 

To summarize, the function f{x) has a local maximum 
at x = 0 and is then strictly monotonically decreasing 
up to Xb ,o (1IV.10I) . Then, for all x > Xb,o the function 
f{x) takes values smaller than or at most equal to the 
values of the envelope function 1 / sin(x), which is strictly 
monotonically decreasing in the considered domain (see 
Figures U and O. Thus, if xi =2 (IIV.7I) is smaller than 
the first boundary point Xb,o (IIV- 10l> . the eigenvalue X 2 
constitutes indeed the second largest eigenvalue. 

Comparing equations (IIV. 71) and (IIV.101) . this is the 
case for TV > 2{k + 1), i.e. for k < TV/2. Numerical 
investigations suggest that the eigenvalue X 2 always con¬ 
stitutes the second largest eigenvalue independent from 
the chosen values for the parameters TV, k and q. How¬ 
ever, monotonicity considerations are not that evident 
for k > N/2. 

The other extremal eigenvalue A mf can not be that eas¬ 
ily assigned to a fixed index. However, arguing similarly 
as for the second largest eigenvalues, it is possible to find 
good estimates for the index at which the smallest eigen¬ 
value 

A mf = min A( nf (IV. 12) 

always occurs. 

The global minimum of f(x) is always located between 
its first two roots x roo t.i and x roo t ,2 in (IIV. 91) . i.e. 


3?root,l 


7T 2tT 

FfT <I - < HT 


% root, 2 • 


(IV.13) 


It thus follows for the index 1. of the smallest eigen¬ 
value: 


TV 

k + 1 


< L - 1 < 


2TV 
k + 1 


(IV.14) 


Therefore, a good estimate for the smallest eigenvalue 
is given by 



topological randomness q 

FIG. 6. Second largest eigenvalues from regular to 
randomized networks. Numerical measurements for undi¬ 
rected (x) and directed (O) networks in comparison with the 
analytical mean field predictions (Eq. (IV.21) . solid lines) as 
a function of q, for different degrees k. The error bars on 
the numerical measurements are smaller than the data points 
(N = 1000, each data point averaged over 100 realizations). 
Adapted from (28j. 


in the following to allow for a consistent analysis for dif¬ 
ferent k. Additionally, we always plot the real part of the 
eigenvalues in the case of directed networks whereas we 
plot A; £ R. in the case of undirected networks. 

As stated in the introduction the principal eigenvalue 
A» (that equals X 2 as shown above) is of special impor¬ 
tance since it dominates the long time dynamics (see e.g. 

0 )- 

For 1 = 2 , eq. fill.2211 simplifies to 


Xf(N,k,q) 




(V.2) 


A mf « Xf 3N +11 , (IV.15) 

L 2(fc + l) +1 l 

where |_£"| denotes the nearest integer to x. 

V. EXTREME EIGENVALUES 

As the offset of each eigenvalue (IIII.22I) equals k , we 
consider the scaled eigenvalues 


In Fig. [6] we compare the typical eigenvalues obtained 
by numerical diagonalization with our analytic prediction 
(IV.21) . The analytic prediction is accurate for both undi¬ 
rected and directed networks, and for all but very small 
relative degrees k/N < 0.05. Moreover, the prediction 
(IV. 21) approximates well the actual dependence of X 2 for 
all but very large q : thus including regular rings, small- 
world and even more substantially randomized network 
topologies. 


A f(N,k,q) 


Ar f (^,M) 

k 


Next we further investigate the scaling behavior of eq. 
(V-l) (IV,21) in dependence on the system’s parameters. 
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FIG. 7. Scaling of the real part of the second largest 
eigenvalue for fixed network size. Numerical measure¬ 
ments for undirected (x) and directed (O) networks in com¬ 
parison with the analytical mean field predictions (Eq. (1 V.2I) , 
solid lines). The analytical approximations obtained from the 
expansions in eq. (1V.3I) are depicted by the dashed lines. The 
error bars on the numerical measurements are smaller than 
the data points, N fixed to 1000. 


A. Approximation for small degrees 



Expanding eq. (IV.21) up to 0(N 2 ) as N —> oo yields 


A? l (N,k,q)~-q- (1 + fc( ^ q))q (V.3) 

(k + 1 )(k + 2) 7t 2 (1 — q) + 6 q{k + 1 — kq) 2 
61V 2 

For <7 = 0 we recover the known approximation for 
symmetric regular ring networks 




(k + l)(fc + 2) 7r 2 
6 N 2 


(V.4) 



The approximation (IV.31) agrees well with eq. (IV.21) up 
to values of k < N/2, but still is a good guide for even 
larger degrees k, cf. Fig. [7) 


B. Scaling with network size 

In order to study the dependence on the network size 
we fix the edge density d = k/N > 0 for large N 1 
which ensures that the network will remain connected. 
This leads to the expression 

a;'V,<,) c -i + ( i (1 _~ ^ (v.5) 


network size N 

FIG. 8. Second largest eigenvalues in dependence on 
edge density and network size, a: Numerical measure¬ 
ments for directed (O) and undirected (x) networks (error 
bars smaller than the data points) in comparison with the 
analytic mean field prediction (IV.51) for N = 2000. b: Asymp¬ 
totic (N —¥ oo) real parts of the second largest eigenvalues A 2 
in dependence on the network size N for fixed edge density 
d = k/N = 0.1 (q -values and symbols as in (a)), c: Asymp¬ 
totic (N —¥ 00 ) real parts of the second largest eigenvalues A 2 
in dependence on the network size N for fixed edge density 
d = k/N = 0.5 (g-values and symbols as in (a)), a and b 
adapted from [28]. 


in the limit N —1 00 . 

Fig- HI confirms the validity of our approximation (IV.51) : 
the second largest eigenvalues A 2 , for both undirected 
and directed networks, in dependence on the edge density 
d for networks of size above about N = 500 nodes are 


predicted well again. Other edge densities than those 
displayed (d = k/N = 0.1 and d = 0.5, Fig. [8] (b),(c)) 
qualitatively yield the same asymptotic behavior. 
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FIG. 9. Smallest eigenvalues from regular to random¬ 
ized networks. Numerical measurements for undirected (x) 
and directed (O) networks in comparison with the analytical 
mean field predictions (Eq. (II V. 121) . solid lines) as a function 
of q, for different degrees k. Dashed lines show the analytic 
estimations of the smallest eigenvalues (Eq. (IIV.15II '). The er¬ 
ror bars on the numerical measurements are smaller than the 
data points (N = 1000, each data point averaged over 100 
realizations). 


C. The smallest eigenvalue 

The smallest eigenvalue A_ defined in Eq. (IIV. 121) is also 
an important indicator for synchronization properties, 
in particular - in combination with the second largest 
eigenvalue - for the synchronizability (see e.g. Hi|). For 
directed networks this refers to the eigenvalue with the 
smallest real part. 

Here, the analytic prediction (IIII.22I) again fits well 
with the actual eigenvalues obtained by numerical diag- 
onalization, cf. Fig. [9] Note also that our estimation 
(IIV.15I) for the smallest eigenvalue agrees well with the 
actual analytic prediction (IIV.12I) . It turns out that the 
analytic prediction is accurate for both undirected and 
directed networks. The prediction (IV. 21) approximates 
well the actual dependence of A_ for small q , thus includ¬ 
ing regular rings and small worlds. The prediction is still 
a good guide for the general dependence of the second 
largest eigenvalue on q , but shows some deviation from 
the numerical results for larger q , i.e. for substantially 
randomized network topologies. 

VI. ANALYTICAL PREDICTIONS VIA 
RANDOM MATRIX THEORY 

To analytically predict the second largest eigenvalues 
for the graph Laplacians of undirected and directed net¬ 
works close to q = 1 (see the shaded area (bottom, right) 
in Fig. 0 we consult Random Matrix Theory [50j (cf. 
also |5ll - l54l |). For a review of synchronization in net¬ 
works with random interactions see e.g. p35l |. 


Firstly, we consider undirected networks associated 
with symmetric matrices. Here, every connection be¬ 
tween a pair of nodes i and j 7 ^ i is present with a given 
probability P. 

Secondly, we consider directed networks associated 
with asymmetric matrices. Here, all nodes have the same 
in-degree = k m . Each of the k m nodes that is con¬ 
nected to node i is independently drawn from the set of 
all other nodes in the network with uniform probability. 

Given a sufficiently large network size N and a suffi¬ 
ciently large k (respectively, a sufficiently large k ln ) we 
numerically find that the set of non-trivial eigenvalues 
resemble disks of radii r' for undirected networks and r 
for directed networks (cf. also BUI). For directed net¬ 
works where the in-degree kl n = k for all nodes i stays 
fixed during the whole rewiring procedure, i.e. all diag¬ 
onal elements are constant, the graph Laplacian is ob¬ 
tained by shifting all eigenvalues of the adjacency matrix 
by —k. For undirected networks there are small devia¬ 
tions from node to node but the average degree equals 
k. However, numerical simulations confirm that shifting 
here again the eigenvalues of the symmetric adjacency 
matrix by the negative average degree —k is feasible. 
Thus, we consider the adjacency matrices in the follow¬ 
ing, A sym for undirected and A asym for directed networks 
and later shift them by —k. 

A. Ensembles of symmetric and asymmetric 

random matrices 

Firstly, consider N x N symmetric matrices A = A T 
with real elements A,j. We constrain the diagonal en¬ 
tries to vanish An = 0 and denote its N eigenvalues by 
A*,. The elements Aij ( i < j) are independent, identi¬ 
cally distributed random variables accor ding to a proba¬ 
bility distribution p(Aij). According to [56h58I | there is 
only one known ensemble with independent identically 
distributed matrix elements that differs from the Gaus¬ 
sian one. Thus there are exactly two universality classes, 
i.e. classes which do not depend on the probability distri¬ 
bution p(Aij), but are determined by matrix symmetry 
only. Every ensemble of matrices within one of these uni¬ 
versality classes exhibits the same distribution of eigen¬ 
values in the limit of large matrices, N —> 00 , but the 
eigenvalue distributions are in general different for the 
two classes. 

The arithmetic mean of the eigenvalues is zero, 

N N 

M i = vE A ‘ = «E' 4 «=° ( vu > 

i=1 i=1 

and the ensemble variance of the matrix elements scales 
like 

- 2 = (4) = £ (VI.2) 

for N 1 and r > 0 being the radii of disks that enclose 
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the set of non-trivial eigenvalues for directed networks 

ii- 

For the Gaussian symmetric ensemble, it is known [50, 
@ that the distribution of eigenvalues Pg™ ss (A) in the 
limit N —^ oo is given by Wigner’s semicircle law 


PgTussW 


^V4r 2 - A 2 if |A| < 2r (VI 3) 
0 otherwise. 


The ensemble of sparse matrices [Sbl . [57l l59l-(62| exhibits a 
different eigenvalue distribution Psp™ se (A) that depends 
on the finite number k of nonzero entries per row and 
approaches the distribution Pg™, ss (A) hi the limit of large 
k such that 


lim PspMse(A) = Pg^ss(A)- (VI.4) 

k — yoo 

It is important to note that in the limit of large N the 
eigenvalue distributions and Pg^ ss depend only on 

the one parameter r, that is derived from the variance of 
the matrix elements (IVI.2I) . 

For real, asymmetric matrices (independent Ajj and 
Aji ), there are no analytical results for the case of sparse 
matrices but only for the case of Gaussian random matri¬ 
ces. The Gaussian asymmetric ensemble yields the dis¬ 
tribution of complex eigenvalues in a disk in the complex 
plane jgH, H3| 

^ otherwise < VL5 > 

where r from Eq. (IVI,2[I is the radius of the disk that is 
centered around the origin. Like in the case of symmetric 
matrices, this distribution also depends only on one pa¬ 
rameter r, that is derived from the variance of the matrix 
elements. 


B. Undirected random networks 


The real symmetric adjacency matrix Al sym is an N x N 
matrix that satisfies and = 0. 

1 J J 1 11 

Furthermore, the matrix elements of A sym are indepen- 
dent up to the symmetry constraint Al^ ym = A??™ ■ They 
are equal to 1 with probability 


= (fc») _ k_ 

IV— 1 ~ N ’ 
and equal to 0 with probability 1 — P. 
Thus, the variance a 2 is given by 

o’-PU-Pj-Ad.*). 


(VI.6) 


(VI-7) 


Therefore, the eigenvalues are located in a disc of ra¬ 
dius 


r' = 2 r 


(VI-8) 


C. Directed random networks 


The real asymmetric adjacency matrix Al asym has ex¬ 
actly k elements equal to one per row. Therefore, its 
elements have a spatial average 

1 N l 

[A 7 m ]==^EA 7 m = ^ ( VL1 °) 

i=i 

and a second moment 

1 N h 

[(-'tf") 2 ] = v E = v - (vi-n) 

j=i 


Thus, the variance 


cr 2 = 




rjasym"|2 _ fZ rZ 

Vi] \ — - jv 2 


(VI.12) 


If we assume that the eigenvalue distribution for di¬ 
rected networks with fixed in-degree is similar to those 
for random matrices [43l . [44| , we obtain a prediction from 
Eq. (IVI.2I) . which yields 



(VI.13) 


for the radius of the disk of eigenvalues centered around 
the origin. 


D. Predictions for the scaled graph Laplacians 

To obtain predictions for the eigenvalues of the appro¬ 
priate graph Laplacian, we have to consider the shift by 
—k (discussed in the beginning of this section) and the 
scaling factor 1/k introduced in eq. (IV.II) . 

Together with eq. (IVI.9I) , the second largest eigenvalues 
for undirected networks close to q = 1 (Fig. [TUI (a)) are 
well predicted by Wigner’s semi-circle law (wsc): 


a r c (v,M) 




(VI.14) 


The real parts of the eigenvalues for directed networks 
close to q = 1 (Fig. [TUI (b)) with eq. (IVI.13I) are with the 
theory of asymmetric random matrices (rmt) 


a r‘WM) 



with 



centered around the origin. 


(VI-9) 



(VI.15) 


Note that A™ sc (iV, k, 1) in eq. (IVI.14I) acquires a posi¬ 
tive value for too small fc-values and a sufficiently large 
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FIG. 10. Analytic prediction of the second largest 
eigenvalues close to q = 1. a: Numerical measurements 
for undirected (x) networks in comparison with the analyti¬ 
cal predictions A™ 11 via Wigner’s semi-circle law (Eq. (IVI.14I) , 
solid lines), for different degrees k. b: Numerical measure¬ 
ments for directed (O) networks in comparison with the an¬ 
alytical predictions A™* from the theory of asymmetric ran¬ 
dom matrices (Eq. (IVI.15I) . solid lines). The error bars on 
the numerical measurements are smaller than the data points 
(N = 1000, each data point averaged over 100 realizations). 
Dashed lines are only a guide to the eye. Taken from [2Sl| . 


network size TV (cf. HI). However, for the A:-values we 
investigated (Fig. [lOl (a) and (b)), the second largest 
eigenvalues are well predicted by both eqs. (1VI.14I) and 

mm - 


VII. SUMMARY AND DISCUSSION 

In this article we have presented and explicated deriva¬ 
tions and extended a simple mean-field rewiring scheme 
suggested recently j28| to derive analytical predictions 
for the spectra of graph Laplacians. The key is re¬ 
placing a stochastic realization of a rewired network at 
given topological randomness by its ensemble-averaged 
network. We achieve this averaging via a two-stage ap¬ 
proach that distinguishes the original outer ring struc¬ 
ture and the originally ’empty’ inner part of a network 
and rewiring probabilities separately. For all q , the re¬ 


sulting average network in particular shares exactly the 
same (average) fraction of links in the original regular 
part of the network as well as in its originally ’empty’ 
part. We derive expressions for the largest nontrivial 
and the smallest eigenvalues, the full spectrum, as well 
as several scaling behaviors. 

We remark that on theoretical grounds, the eigenvalue 
spectrum of the resulting average network in general is 
not equal to the average of the spectra of the individual 
stochastic network realizations. Yet, systematic numer¬ 
ical checks confirm that the mean field approximation 
introduced is accurate as long as q is sufficiently below 
one. In the limit q —> 1, we derive the eigenvalue spectra 
based on random matrix theory which become exact in 
the limit of infinitely large networks, TV —> oo. 

Although the mean field rewiring is undirected, eigen¬ 
values for directed networks are approximated more ac¬ 
curately and in a wider range of q- values, which is in 
particular related to the fact that the predictions for the 
undirected second-largest eigenvalues at q = 1 are larger 
in real part than the directed ones, while all the mean 
field eigenvalues converge to —1 at q = 1. For ‘small’ de¬ 
values the mean field approximation becomes less accu¬ 
rate, which may be due to the fact that the ring structure 
is destroyed more easily while rewiring. Additionally, the 
bulk spectra spread much more drastically with q than 
for larger k- values. 

Furthermore, note that the analysis of the mean field 
spectrum presented here can principally not be extended 
to the Laplacian eigenvectors as these are independent 
of the mean field Laplacian’s elements Cj IIII. 81 just be¬ 
cause of the circulant structure of the mean field graph 
Laplacian. Consequently, the eigenvectors are the same 
and always non-localized, independent of the system’s 
parameters TV, k and q. Studying distributed patterns of 
relaxation processes and potential localization phenom¬ 
ena thus requires access to eigenvectors beyond the mean 
field approximation. Studies of the Laplacian eigenvec¬ 
tors are rare although there are fascinating results as 
well. For instance, the discrete analogs of solutions of 
the Schrodinger equation on manifolds can be investi¬ 
gated on graphs (cf. e.g. 0). 

The simple mean field approach presented above still 
substantially reduces computational efforts when study¬ 
ing randomized (regular or small-world) network models. 

Generalizing our mean field approach to higher di¬ 
mensions and/or to other rewiring approaches, as for in¬ 
stance, relevant for neural network modeling [66|, it will 
serve as a powerful tool to gain new insights into the re¬ 
lations between structural and dynamical properties of 
complex networks. 
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